﻿//******************************************************************
//  caja.h
//  Objeto 'caja de simulacion'  para MDDipFlow.cpp
//	Grupo de Fisica de Fluidos y Biocoloides - Univers. de Granada
//  ©2011, Dr. Juan Carlos Fernandez Toledano
//******************************************************************
class caja
{
private:

public:
    ///Lado de la caja
	double lbox[3], lboxh[3];
	///Array de particulas
	particula *part;
	///Numero de particulas
	int npart;
	///Concentracion en volumen
	double phi;
	///Concentracion en masa
	double rho;
	///Temperatura de caja
	double temp;
	///Temperatura de referencia
	double temp0;
	///Volumen
	double vol;
	///Distancia de corte de interaccion (NO UTILIZADO)
	double rcut;
	///Distancia a pares entre particulas
	double dr[3];
	///Distancia a pares entre particulas sin condiciones periodicas
	double drnb[3];
	///Presion
	double presion;
    ///Altura de la caja
    double h0;
	///Energia potencial
	double epoten;
	///Virial
	double virial;
	///Energia cinetica
	double ekin;
	///Control de la velocidad (termostato)
	double alpha;
	///Paso de la simulacion
	int nacc;
	///Velocidad media de las particulas de la caja
	double vel[3];
	///Elementos del 'microscopic stress tensor' asociado a la componente cinetica
	double svxy, svxz, svyz, svzz;
	///Elementos del 'microscopic stress tensor' asociado al potencial de interaccion
	double spxy, spxz, spyz, spzz;
	///Campo magnetico
	double u0, u30;
	///Parametros de la simulacion
	parametros param;
	///shear rate
	double gamma_punto;
	///Velocidad de la pared
	double vel_shear;

	caja(parametros par)
	{
		///Empleamos paso de tiempo Dt para generar las posiciones posOld
		///de las particulas empleando las velocidades
		param=par;
		if(param.rest==0)
		{
			///Iniciamos la simulacion desde el principio
			u0=param.u0;
			u30=3.0*u0;
			h0=param.h0;
			npart=param.npart;
			//gamma_punto=param.shear;
			vel_shear=gamma_punto*h0;
			temp0 =1.0;
			phi=param.phi;
			rcut = param.rcut;
			rho = phi*6.0/pi;
			vol=pi*npart/phi/6.;
			//vol=pi*npart/phi/6.0;
			#ifdef CONFINED
			if (h0>0.){
                ///La altura de la caja es tamayo caja menos 2a
                lbox[0]=pow(vol/h0,1./3.);
                lbox[2]=h0*lbox[0];
                lbox[1]=lbox[0];
			}
			else{
			    lbox[0]=pow(vol,1./3.);
			    lbox[1]=lbox[2]=lbox[0];
			}
			#endif
			#ifndef CONFINED
			lbox[0]=pow(vol,1./3.);
            lbox[1]=lbox[2]=lbox[0];
            #endif
            for(int i=0;i<3;i++) lboxh[i]=0.5*lbox[i];
			part = new particula[npart];
			SetRandomPart();
			initForces();
			initVelocities(param.Dt);

		}
		else
		{
			///Reiniciamos una simulacion desde archivo existente
			loadConfig();
			for(int i=0;i<3;i++) lboxh[i]=0.5*lbox[i];
		}
	}

	~caja()
	{
		delete [] part;
	}

	///Imprime la configuracion de reinicio
	void printConfigBackup();
	///Coloca las particulas de forma aleatoria en una caja
	void SetRandomPart(void);
	///Inicializa las fuerzas
	void initForces(void);
	///Inicializa las velocidades y calcula posOld
	void initVelocities(double Dt);
	///Aplica condiciones contorno periodicas a particula ind
	void boundary(int ind);
	///Calcula la distancia cuadratica entre particulas
	inline double distancia2(int ind1,int ind2);
	///Calcula la distancia entre particulas
	inline double distancia(int ind1,int ind2);
	///Imprime el archivo XYZ de la configuracion de particulas
	void printXYZfile(char namefile[30]);
	///Se imprime el tensor de esfuerzos microscopicos
	void printShearStressTensor();

	///Carga de una configuracion para el reinicio de la simulacion
	void loadConfig(void);
	///Importa diametros
	double importDiametros(double *diam, int npart);
	void expand();
	void contract();
};

void caja::loadConfig(void)
{
	ifstream MiArchivo ("system_ini.dat",ios::in);
	char line[150];
	double dummy2;
	double rr;
	if(MiArchivo.fail())
	{
		cout<<"No se pudo encontrar el archivo de reinicio "<<endl;
		exit(1);
	}
	if (MiArchivo.is_open())
	{
        MiArchivo.getline(line, sizeof(line));
        MiArchivo.getline(line, sizeof(line));
        MiArchivo.getline(line, sizeof(line));
        MiArchivo.getline(line, sizeof(line));
        MiArchivo>>npart>>phi>>rcut>>dummy2>>dummy2>>dummy2>>dummy2>>lbox[0]>>lbox[1]>>lbox[2]>>dummy2>>dummy2;
		///Construimos las cosas
		u0=param.u0;
		u30=3.0*u0;
		npart=param.npart;
		cout<<npart<<endl;
		temp0 =1.0;
		phi=param.phi;
		rcut = param.rcut;
		rho = phi*6.0/pi;
		part = new particula[npart];
		initForces();
		vol=lbox[0]*lbox[1]*lbox[2];
		lboxh[0]=lbox[0];
		lboxh[1]=lbox[1];
		lboxh[2]=lbox[2];
		///Fin de la construccion

        MiArchivo.getline(line, sizeof(line));
		MiArchivo.getline(line, sizeof(line));
		MiArchivo.getline(line, sizeof(line));
		MiArchivo.getline(line, sizeof(line));
		MiArchivo.getline(line, sizeof(line));
		MiArchivo.getline(line, sizeof(line));
		///Importamos las posiciones y velocidades
		for(int i=0;i<npart;i++) {
			for(int j=0;j<3;j++){
				MiArchivo>>rr;
				part[i].pos[j]=rr;
			}
			for(int j=0;j<3;j++){
			    MiArchivo>>rr;
				part[i].vel[j]=rr;
			}
//			for(int j=0;j<3;j++){
//			    MiArchivo>>rr;
//				part[i].force[j]=rr;
//
//			}
			boundary(i);
			#ifdef POLYDISPERSE
                MiArchivo>>rr;
                part[i].radio=rr;
			#endif
			part[i].mass=pow(part[i].diametro,3.0);
		}
		MiArchivo.close();
	}
}

void caja::printConfigBackup()
{
	char name[50];
	sprintf(name, "%s","system_bk.dat");
	ofstream MiArchivo (name,ios::out);
	//MiArchivo.precision(4);

	if (MiArchivo.is_open()) {
		MiArchivo.setf(ios::fixed);
		MiArchivo<<"Parametros de la simulacion"<<endl;
		MiArchivo.fill('*');
		MiArchivo << setw(122) << ""<<endl;
		MiArchivo.fill(' ');
		MiArchivo<<"Npart"<<setw(10)<<"phi"<<setw(10)<<"rcut"<<setw(10)<<"Dt"<<setw(10)<<"restart"
			<<setw(10)<<"trelaj1"<<setw(10)<<"bx"<<setw(10)<<"by"<<setw(10)<<"bz"<<setw(10)<<"tiempo"<<setw(10)<<"paso"<<endl;
		MiArchivo.fill('-');
		MiArchivo << setw(122) << ""<<endl;
		MiArchivo.fill(' ');
		MiArchivo<<npart<<setw(10)<<phi<<setw(10)<<rcut<<setw(10)<<param.Dt<<setw(10)<<param.rest
			<<setw(10)<<lbox[0]<<setw(10)<<lbox[1]<<setw(10)<<lbox[2]<<setw(10)<<param.Dt*nacc<<setw(10)<<nacc<<setw(10)<<endl;
	}
	MiArchivo<<endl;
	MiArchivo<<"Posiciones y velocidades"<<endl;
	MiArchivo.fill('*');
	MiArchivo << setw(122) << ""<<endl;
	MiArchivo.fill(' ');
	MiArchivo<<"rx"<<setw(20)<<"ry"<<setw(20)<<"rz"<<setw(20)<<"vx"<<setw(20)<<"vy"<<setw(20)<<"vz"<<endl;
	MiArchivo.fill('-');
	MiArchivo << setw(122) << ""<<endl;
	MiArchivo.fill(' ');
	for(int i=0;i<npart;i++) {
	    MiArchivo<<scientific<<part[i].pos[0]<<setw(20)<<part[i].pos[1]<<setw(20)<<part[i].pos[2]<<setw(20)<<part[i].vel[0]<<setw(20)<<part[i].vel[1]<<setw(20)<<part[i].vel[2];
	    //MiArchivo<<setw(20)<<part[i].force[0]<<setw(20)<<part[i].force[1]<<setw(20)<<part[i].force[2]<<endl;
	    #ifdef POLYDISPERSE
            MiArchivo<<setw(20)<<part[i].radio;
	    #endif
        MiArchivo<<endl;
	}
}


void caja::SetRandomPart(void)
{
    int nok=0;
    bool good;
    double r2;

    #ifdef POLYDISPERSE
        double poly=0.3;
        double Vpart=0.0;
        double u1,u2;
        ///Calculamos el tamanyo de caja para la polidispersidad considerada
        //Vpart = importDiametros(diametro, npart);
        for(int i=0;i< npart;i++)
        {
            u1=double(rand())/RAND_MAX;
            u2=double(rand())/RAND_MAX;
            //Aqui introduciremos la distribucion de tamaños
            part[i].radio=exp(poly*sqrt(-2.*log(u1))*cos(2.*pi*u2));
            Vpart+=pi*pow(part[i].radio,3.0)/6.0;
            part[i].mass=pow(part[i].radio,3.0);
        }
        vol=Vpart/phi;
        #ifdef CONFINED
            if (h0>0.){
                lbox[0]=pow(vol/h0,1./3.);
                lbox[2]=h0*lbox[0];
                //lbox[0]=sqrt(vol/lbox[2]);
                lbox[1]=lbox[0];
            }
            else{
                lbox[0]=pow(vol,1./3.);
                lbox[1]=lbox[2]=lbox[0];
            }
        #endif
        #ifndef CONFINED
            lbox[0]=pow(vol,1./3.);
            lbox[1]=lbox[2]=lbox[0];
        #endif
    #endif

    #ifndef POLYDISPERSE
        double Vpart=0.0;
        for(int i=0;i< npart;i++)
        {
            part[i].radio=1.0;
            Vpart+=pi*pow(part[i].radio,3.0)/6.0;
            part[i].mass=pow(part[i].radio,3.0);
        }
        vol=Vpart/phi;
        #ifdef CONFINED
            if (h0>0.){
                lbox[0]=pow(vol/h0,1./3.);
                lbox[2]=h0*lbox[0];
                //lbox[0]=sqrt(vol/lbox[2]);
                lbox[1]=lbox[0];
            }
            else{
                lbox[0]=pow(vol,1./3.);
                lbox[1]=lbox[2]=lbox[0];
            }
        #endif
        #ifndef CONFINED
            lbox[0]=pow(vol,1./3.);
            lbox[1]=lbox[2]=lbox[0];
        #endif
    #endif

	for(int i=0;i< npart;i++)
	{
	    part[i].pos[0]=lbox[0]*(double(rand())/RAND_MAX);
        part[i].pos[1]=lbox[1]*(double(rand())/RAND_MAX);
        part[i].pos[2]=(lbox[2]-2.*part[i].radio)*(double(rand())/RAND_MAX)+part[i].radio;

		boundary(i);
		///Comprobamos solapamiento con otras particulas
		good=true;
		for(int j=0;j<nok;j++){
		    ///Recorremos todas las particulas que ya esten metidas en la lista
		    if(good){
		        r2=distancia2(i,j);
		        #ifdef POLYDISPERSE
                if(sqrt(r2)<0.5*(part[i].radio+part[j].radio)){
                    ///Hay solapamiento entre las particulas i y j. No aceptamos la particula i
                    good=false;
                    i--;
                }
                #endif

                #ifndef POLYDISPERSE
                if(r2<1.2){
                    ///Hay solapamiento entre las particulas i y j. No aceptamos la particula i
                    good=false;
                    i--;
                }
                #endif
		    }
		}
		if(good){
            ///Si pasa por aqui la particula no se solapa con ninguna otra --> se acepta
            nok++;
		}
	}
	cout<<"Particulas colocadas en la caja: "<<endl;

}

void caja::initForces(void)
{
	for (int i=0;i<npart;i++) {
		for (int j=0;j<3;j++){
			part[i].force[j]=0;
		}
	}
}

void caja::initVelocities(double Dt)
{
	//double v0=sqrt(3.0*kB*temp);
	double v02=0.0;
	double dv[3];
	double enerkin,enerkin2;

	for (int i=0;i<npart;i++) {
	    for (int j=0;j<3;j++) {
			part[i].vel[j]=0.0;
		}
	}
}

void caja::expand()
{
    for(int i=0;i<npart;i++)
    {
        for(int j=0;j<3;j++)
        {
            part[i].pos[j]=(part[i].pos[j]+0.5*lbox[j])*0.5;
        }
    }
}

void caja::contract()
{
    for(int i=0;i<npart;i++)
    {
        for(int j=0;j<3;j++)
        {
            part[i].pos[j]=2.*part[i].pos[j]-0.5*lbox[j];
        }
    }
}

void caja::boundary(int ind)
{
//    int ndim=3;
//    #ifdef CONFINED
//    ndim=2;
//    #endif
//	for(int i=0;i<ndim;i++)
//	{
//		part[ind].pos[i]-=lbox[i]*anint(part[ind].pos[i]/lbox[i]);
//	}
    //if(part[ind].pos[0]<0)//lbox[0])
    //    cout<<dround(part[ind].pos[0]/lbox[0])<<"  "<<part[ind].pos[0]<<"  "<<lbox[0]<<endl;
	if(part[ind].pos[0]>lbox[0]) part[ind].pos[0]-=lbox[0];
	if(part[ind].pos[0]<0.) part[ind].pos[0]+=lbox[0];

	if(part[ind].pos[1]>lbox[1]) part[ind].pos[1]-=lbox[1];
	if(part[ind].pos[1]<0.) part[ind].pos[1]+=lbox[1];


	#ifndef CONFINED
	if(part[ind].pos[2]>lbox[2]) part[ind].pos[2]-=lbox[2];
	if(part[ind].pos[2]<0.) part[ind].pos[2]+=lbox[2];
	#endif

}

inline double caja::distancia2(int ind1,int ind2)
{
	///Calcula la distancia entre las particulas i y j
	///Consideramos condiciones de contorno periodicas en el plano xy
	///pero no en el z donde se encuentran las paredes

	double dis2=0;
//    dr[0]=part[ind2].pos[0]-part[ind1].pos[0];
//    dr[0]-=lbox[0]*anint(dr[0]/lbox[0]); //Entero + proximo manteniendo signo
//    /*Y*/
//    dr[1]=part[ind2].pos[1]-part[ind1].pos[1];
//    dr[1]-=lbox[1]*anint(dr[1]/lbox[1]); //Entero + proximo manteniendo signo
//    /*Z*/
//    dr[2]=part[ind2].pos[2]-part[ind1].pos[2];
//    #ifndef CONFINED
//    dr[2]-=lbox[2]*anint(dr[2]/lbox[2]); //Entero + proximo manteniendo signo
//    #endif
//
//	dis2=dr[0]*dr[0]+dr[1]*dr[1]+dr[2]*dr[2];
//cout<<lboxh[0]<<"  "<<lboxh[1]<<"  "<<lboxh[2]<<endl;
	/*X*/
	dr[0]=part[ind2].pos[0]-part[ind1].pos[0];
	//if(abs(dround(dr[0]/lbox[0]))==1) cout<<round(dr[0]/lbox[0])<<"  "<<dr[0]<<"  "<<lbox[0]*0.5<<"  "<<dr[0]-dround(dr[0]/lbox[0])*lbox[0]<<endl;
    //dr[0]-=round(dr[0]/lbox[0])*lbox[0];
    if(dr[0]>lboxh[0]) dr[0]-=lbox[0];
	if(dr[0]<-lboxh[0]) dr[0]+=lbox[0];
	//dr[0]-=lbox[0]*anint(dr[0]/(lbox[0])); //Entero + proximo manteniendo signo
    /*Y*/
    dr[1]=part[ind2].pos[1]-part[ind1].pos[1];
    //dr[1]-=round(dr[1]/lbox[1])*lbox[1];
    //dr[1]-=lbox[1]*anint(dr[1]/(lbox[1])); //Entero + proximo manteniendo signo
    if(dr[1]>lboxh[1]) dr[1]-=lbox[1];
	if(dr[1]<-lboxh[1]) dr[1]+=lbox[1];
    /*Z*/
    dr[2]=part[ind2].pos[2]-part[ind1].pos[2];
    #ifndef CONFINED
    //dr[2]-=round(dr[2]/lbox[2])*lbox[2];
    if(dr[2]>lboxh[2]) dr[2]-=lbox[2];
    if(dr[2]<-lboxh[2]) dr[2]+=lbox[2];
    #endif

	dis2=dr[0]*dr[0]+dr[1]*dr[1]+dr[2]*dr[2];

	return	 dis2;
}

inline double caja::distancia(int ind1,int ind2)
{
	return	 sqrt(distancia2( ind1,ind2));
}

void caja::printXYZfile(char namefile[30])
{
	char name[50];
	sprintf(name, "%s", "representacion.dat");
	ofstream MiArchivo (name,ios::out|ios::app);
	MiArchivo.precision(5);
	if (MiArchivo.is_open())
	{
		MiArchivo << npart<<endl;
		MiArchivo << lbox[0]<<"  "<<lbox[1]<<"  "<<lbox[2]<<"  "<<nacc<<endl;
		for(int i=0;i<npart;i++)
		{
			MiArchivo<<part[i].pos[0]<<"  "<<setw(8)<<part[i].pos[1]<<"  "<<setw(8)<<part[i].pos[2]<<"  "<<setw(8)<<part[i].radio<<endl;
		}
		MiArchivo.close();
	}
}

void caja::printShearStressTensor()
{
	char name[50];
	sprintf(name, "%s", "stressTensor.dat");
	ofstream MiArchivo (name,ios::out|ios::app);
	if (MiArchivo.is_open())
	{
		MiArchivo << nacc*param.Dt<<"  "<<svxy+spxy<<"  "<<svxz+spxz<<"  "<<svyz+spyz<<"  "<<svzz+spzz<<endl;
		MiArchivo.close();
	}
}

